#include "PO_MoM.h"
#include "FillMatrix.h"
#include <iostream>
using namespace blitz;

void FillV()
{
	V.resize(iWire);
    V(iWire/2) = 1;
}
void FillZ()
{
	Array<dComplex,2> L_J_E(iWire,iEdge), L_I_E(iWire,iWire);
	Z.resize(iWire,iWire);
	Tau.resize(iEdge,iWire);
	for (int n=0; n<iWire; n++)
		for (int k=0; k<iEdge; k++)
		{
			Tau(k,n) = fun_Tau(k,n);
			L_J_E(n,k) = fun_L_J_E(n,k);
		}
	for (int m=0; m<iWire; m++)
		for(int n=m; n<iWire; n++)
		{
			L_I_E(m,n) = fun_L_I_E(m,n);
			L_I_E(n,m) = L_I_E(m,n);
		}

	firstIndex i;
	secondIndex j;
	thirdIndex k;
	Z = L_I_E - sum(L_J_E(i,k) * Tau(k,j),k) * DeltaLength;
}
dComplex fun_Tau(int k, int n)
{
	Position3D Norm = GetEdgeNormal(k);
	Position3D R = GetR(k,n);
	if (dot(R,Norm) < 0)
	{
		Complex3D tmp1 = cross(Nabla_G(R),Complex3D(UnitVector))*DeltaLength/(4*PI);
		Complex3D tmp2 = cross(Complex3D(Norm),tmp1);
		return dot(tmp2,T(k));
	}
	else return 0;
}
dComplex fun_L_I_E(int m, int n)
{
	return dComplex(0,1)*Omega*Mu/(4*PI)*DeltaLength*DeltaLength*G(vWire(m)-vWire(n))
		-dComplex(0,1)/(Omega*Epsilon*4*PI)*(G(vWire(m)-vWire(n))-G(vWire(m)-vWire(n)+DeltaLength*UnitVector)-G(vWire(m)-vWire(n)-DeltaLength*UnitVector)+G(vWire(m)-vWire(n)));
// 	return dComplex(0,1)*Omega*Mu/(4*PI)*DeltaLength*DeltaLength*G(vWire(m)-vWire(n))
// 		+dComplex(0,1)/(Omega*Epsilon*4*PI)*dot(Position3D(0,0,DeltaLength),(Nabla_G(vWire(m)-vWire(n)+Position3D(0,0,DeltaLength/2))-Nabla_G(vWire(m)-vWire(n)-Position3D(0,0,DeltaLength/2))));
}
dComplex fun_L_J_E(int m, int k)
{
	Complex3D tmp1 = (Nabla_G(vWire(m)-vPatchCtr(Edge(k,2)))-Nabla_G(vWire(m)-vPatchCtr(Edge(k,3))))*dComplex(0,-Edge_length(k)/(4*PI*Epsilon*Omega));
	Complex3D tmp2 = (GetRho(k,0)*G(vWire(m)-vPatchCtr(Edge(k,2)))+GetRho(k,1)*G(vWire(m)-vPatchCtr(Edge(k,3))))*dComplex(0,-Edge_length(k)*Omega*Mu/(8*PI));
	return dot(tmp1+tmp2,UnitVector);
}
dComplex G(Position3D R)
{
	double RLength = sqrt(dot(R,R));
	if (RLength < 1e-7) return SingularG();
	else return dComplex(cos(-WaveNumber*RLength),sin(-WaveNumber*RLength))/RLength;
}
Complex3D Nabla_G(Position3D R)
{
	double RLength = sqrt(dot(R,R));
	dComplex tmp = dComplex(cos(-WaveNumber*RLength),sin(-WaveNumber*RLength))/(RLength*RLength*RLength)*dComplex(-1,-WaveNumber*RLength);
	return R*tmp;
}
Position3D GetEdgeNormal(int k)
{
	Position3D tmp = (vPatchNorm(Edge(k,2))+vPatchNorm(Edge(k,3)))/2;
	return tmp/sqrt(dot(tmp,tmp));
}
Position3D GetR(int k,int n)
{
	Position3D tmp=(vNode(Edge(k,0))+vNode(Edge(k,1)))/2;
	return tmp-vWire(n);
}
Position3D GetRho(int k, int sign)
{
	if (sign == 0)
	{
		int pat = Edge(k,2);
		int node = Patch(pat,0)+Patch(pat,1)+Patch(pat,2)-Edge(k,0)-Edge(k,1);
		return vPatchCtr(pat)-vNode(node);
	}
	else
	{
		int pat = Edge(k,3);
		int node = Patch(pat,0)+Patch(pat,1)+Patch(pat,2)-Edge(k,0)-Edge(k,1);
		return vNode(node)-vPatchCtr(pat);
	}
}
dComplex SingularG()
{
	return dComplex(2/DeltaLength*log(DeltaLength/WireRadius),-WaveNumber);
}
